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1. Introduction 



Let II = [^j+fc]^fc=:o be a be a Hankel operator defined by its rational symbol p = tt/x, where 

n— 1 n 



j=0 J=0 



( 1 . 1 ) 



We assume that tt and x have no common zeros. The elements r]j of H are then given by 



( 1 . 2 ) 



In order to simplify our presentation, we assume that the zeros {^k)k=i X distinct. How our 
formulas need to be modified in order to remove this assumption is discussed in Remark 1.1 below. 
Hence p ha^ a partiaJ fraction decomposition 






^=1 



(1.3) 



Expansion of the right hand side of (1.3) into a geometric series, and comparison with (1.2), yields 

n 

= (1.4) 

A:=l 

We now express (1.4) in matrix form. Let 

A : = diag[ai,a2,... a„] e 

A:= d>aj^[Ai,A2,... A„] c 

and introduce the Vandermonde matrix 

:= [V,]°lo f 

V,- := VoA^'*, j>l. 

H = VA V^. 



Define 

where 

Then (1.4) can be written eis 



(1.5) 

( 1 . 6 ) 

(1.7) 

( 1 . 8 ) 
(1.9) 

( 1 . 10 ) 



Let denote the vector space <C°° equipped with the Euclidean norm. 
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Proposition 1.1. bounded jA^I < 1 for 1 < A: < u. 

Proof. The proposition holds independent of the multiplicity of the A^. In the present proof we 
assume that the A^- are distinct. The proof for confluent is commented on in Remark 1.1. 

Let €i = be the axis vector wdth cq = \ . Then 

h = •= -A/ei G r/j — > 0 as j — > oo => 

|A^| < 1 for 1 < A: < n, 
where the last implication follows from (14). 

Conversely, assume that \X^\ < 1 for 1 < A: < n. Then by (1.8) - (1.10) w^e obtain 

oo 

ll/^lb < ll^lbllV^II^ < ||>4||2||^o||2"|lE(A''^ni2 

= II^I|2||v^oII^II(/-(a^a)-)'^||^ . 

We assume henceforth that \X^\ < 1 for 1 < A: < n. Introduce 

U:^VVo\ (1.11) 

Ho:=VoAV^^. (1.12) 

Then Hq has rank n. We note , by comparing (1.12) with (l.lO), that Hq is the leading principal 
n X n submatrix of H. From (l lO) - (112) it follows that 

H = UHoU'^. (1.13) 

The leading n X n submatrix of U is 7^, the n x n identity matrix. U therefore is of rank n and 

can be factored 

U = QR, 

where Q^Q = In ^.nd 7? is a nonsingular right triangular matrix. We obtain 

a+{ll) = o^{QRUolfQ'^) = o{RlhR^), (1.14) 

w'here o denotes the set of singular values and denotes the subset of the positive ones. 

The n X n matrix RIIqR^ is complex symmetric. Takagi [Tal], [Ta2] showed the existence of a 
complex symmetric singular value decomposition 

RI/oR'^ = \V( i: = dt'aglai,a2, ... (1.15) 

where W^W = /„ and <7j > 0 arc the singular values of RH^R^ . In Section 2 we present an 
elementary proof of the existence of this decomposition. Let IV = [w],ui 2 ,... »f„], Wj c (T". Then 
(1.15) can be written as the Takagi singular value problem 

RHqR^Wj = WjOj, = 8jk, 1 < J, (I 16) 
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where the bar denotes complex conjugation and is Kronecker’s 8 function. The problems (1.15) 
- (1-16) could be solved by the algorithm described in [BGG], but this would require RHqR^ to be 
explicitly computed. In order to avoid these matrix multiplications we let vy := R^ xVj and obtain 
from (1.16) the generalized Takagi singular value problem 

HqVj = (R^ R)~^Vj<7j, {R^ R)~^Vk = Sjk, 1 < j,k < n. (1-17) 

The solution of (1-17) requires (R^ R)~^ to be known. In Section 3 we show that 

{R’^R)-^ = I - BoB^, (1.18) 

where Bq € is a triangular Toeplitz matrix. The elements of Bq and Hq can be determined 

from the coefficients of tt and x i*^ 0{n log n) arithmetic operations by the fast Fourier transform 
(FFT) method. This is demonstrated in Section 4. Section 5 shows that 

= Ti,Mo e (1.19) 

where Tj and Mq are Toeplitz matrices, and describes a numerical scheme for the computation 
of this factorization from (1.16) in O(n^) arithmetic operations. We also present a Hermitian 
factorization of R^ R into n X n triangular matrices. 

The f 2 u:torization (1.19) may be of interest for the numerical solution of (1.17). Assume that 
the coefficients of and x ^re real valued. Then //q, {R^ R)~^ e and (1.17) reduces to 

a generalized symmetric eigenvalue problem. The Lanczos method ([Pa, Section 15.11], [ER]) 
would appear suitable for solving this eigenproblem for the following reason. Let C e be a 

Hankel or Toeplitz matrix and let v € be arbitrary. It is well known that Cv can be computed 
in 0(n log n) arithmetic operations using FFTs. Hence Hqv^ {R^ R)~’^v and {R^ R)v can be 
computed in 0(n log n) arithmetic operations, where we use (1.18) - (1.19). Each iteration of the 
Lanczos algorithm given in [Pa, p.324| therefore requires only 0(n log n) arithmetic operations. 

The computation of singular values of H is important in Hankel norm approximation problems of 
systems theory, such as the model reduction problem [Glj. The approximation of functions by the 
Caratheodory - Fejer method yields another application [Gu], [Trj. 

Other methods for reducing the singular value problem for 77 to a finite dimensional one have been 
described by Rung and Gutknecht [Gu] and Young [Yo]. These methods, however, do not preserve 
symmetry. Moreover, Young’s approach requires generally O(n^) arithmetic operations to compute 
the matrices required. 



Remark 1.1. Formulas (1.3) - (1.8) and the proof of Proposition 1.1 require 
restriction can be removed. Assume first that Ai = A 2 = ... = A^. Then (1.3) 
replaced by 



0(A) =: x: 



^ = 1 



{X-Xif' 



distinct This 
- ( 1 . 4 ) have to be 



(1.3') 






n CX3 . 



k=l 



3-0 



(1.4') 



In ( 1 . 5 ) A has to be substituted by the upper triangular Hankel matrix 

= K + f ap:=0, p>n. 
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The matrix A in (1.6) has to be replaced by the Jordan matrix with all diagonal elements equal to 
Ai and all superdiagonal elements equal to one. The matrix Vq in (1.7) need be replaced by the 
confluent Vandermonde matrix. For instance, we obtain for n = 3 





Ai 1 




1 




A = 


Aj 1 


Vo = 


Ai 


1 




'^1 ^ 




A^ 

L^i 


— 1 



With Ay A and Vq modified as described, we define Vj and V by (1.8) - (1.9), U by (1.11) and 
Ho by (1.12). Then (1.10) and (1.13) hold and Hq is the leading principal n x n submatrix of H. 
Also (1.14) - (1.19) remain valid. Proposition 1.1 can be shown by replacing (1.4) by (1.4'), and 
by bounding the sum 

oo 

j=0 

w^here A now is a Jordan matrix. This sum is bounded if \X\ | < 1, and the proposition remains 
valid. 

In general, when the A/^ are of arbitrary multiplicity, A in (1.5) hats to be replaced by a block 
diagonal matrix, where each block is an upper triangular Hankel matrix. The blocks are of the 
same sizes as the multiplicities of the A^c , and the number of blocks equals the number of distinct A^ . 
A in (1.6) is replau:ed by a Jordan matrix with Jordan boxes of the same sizes as the multiplicities 
of the A/c) and the number of boxes equal to the number of distinct A^. Vo (l-*^) is replaced by 
an appropriate confluent Vandermonde matrix. With these changes (1.10) - (1.19) are valid, and 
so is Proposition 1.1. We omit the details since the numerical computations are independent of the 
multiplicity of the A^^- • 
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2. The Symmetric Singular Value Decomposition 



In this section we present an elementary proof of Takagi’s theorem, i.e. we 
a symmetric singular value decomposition of a complex symmetric matrix, 
and define A, B e by C := yl + iB, i := Then A = A^and B = 



show the existence of 
Let C = C^ € 

B^ , so the matrix 



C := 



A 

B 



B' 

-A 



is reed and symmetric. Let be the positive eigenvalues of C and form 



Let 



with 



and 



E := diag{ai,a2,...,ar]. 



A 


b' 




'u' 




'u' 


B 


-A 




V 







U,Ve 

U'^U + V'^V = Ir. 



Write (2.1) as 

<AU + BV = UT, 
\BU - AV = VE 

and note that (2.2) also can be written as 



( 2 . 1 ) 



( 2 . 2 ) 



i.e. 



with 



AV 


+ B{-U) = 


V{-E) 


BV 


- A{-U) = 


(-t/)(-E), 


A 


Bl [ Fl 




B 


-A -U 



V^V + (-[/f(-U) = Ir. 



Hence C has at least r negative eigenvalues. We could also have let Cj be the negative eigenvalues of 
C and then (2.3) would have given us positive ones. We therefore may assume that ±ai , ±^ 2 , ±<7^ 
are all the nonzero eigenvalues of C. 



Since eigenvectors associated with distinct eigenvalues of a 
we have 



0 = 



u 

V 



v'^u 



real symmetric matrix are orthogonal, 
- U'^V. 



The spectral resolution of C is thus 



A B' 




'U V' 




E 




■[/^ 1/^' 


B -A 




y -U 




-E 
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which yields 

( A = [/Sf7^ - vi:v^ 

[ D = VEU'^ + UEV^ 



Therefore 

where 


C = A + iB = UEU"^ - VEV'^ + i[VEU'^ + UEV^) 

r 

= [U -f t'V)E(U^ + .V^) = WEW'^ = Y^OkWwl 

fc=l 


Moreover 


U + iV =: W = [u;i,it;2, Wr], Wfcf ■ 



= ([/^ - iV'^){U + iV) = {u'^u + V'^V) + i{u'^V - V'^U) 

If r < n then one may replace E by 





Eo ;= diag[oi,a'2,...,(Tr,0,...,0] ( 


and W by 


IV'o := [ti;i , W2, ...,u;r,u;r + i, f 


where u;r-fi, 


C" are chosen so that W^Wq — 1,^. ■ 






3. A Simple Expression for {R^ R) ^ 

In this section we derive (1.18). Introduce the Frobenius matrix 



F:= [e2,e3,...,e„,-/] 

where 

ey : = eR", 2 < j < n, (3.1) 

/ : = [xo,Xi,-.Xn-i]^f 

Then F is the companion matrix of x ^.nd 



F^Vo = VoA. (3.2) 

Throughout this section Vq and A are defined by (16) - (1.7) if the are distinct. For confluent 
Xk we modify Vq and A according to Remark 1.1. The following lemma shows that 

G := R^R (3.3) 

satisfies a Stein equation. This will enable us to obtain a simple expression for by an application 
of the Sherman-Morrison- Woodbury formula. 

Lemma 3.1. G is the unique solution of the Stein equation 

X - = Iny Xe (3.4) 

Proof. By (1.8), (1.9) and (1 11) we obtain 

cx> 

R^R = U^U = ^ Vo-"(A"^)"Vo"VoA"^Vo“\ (3.5) 

k=0 



and (3.2) yields now 

CO 

G = 

fc = 0 



(3.6) 



The series in (3.5) - (3.6) converge because \Xk\ < 1 for all k. Substitution of (3.6) into (3.4) 
shows that G solves (3.4). The unicity follows from \Xk\ < 1 for all k. The latter can be seen by a 
similarity transform of F^ to Schur triangular form. « 

Introduce the cyclic downshift operator in 

E := ( 62 , 63 ...., e„,ei] 6 



where 



Let 



6j := f R^"- 

i ~ (xo,Xi.-.X«,0,0,...,0]^€ 



(3.7) 
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and define the Toeplitz matrix T of parallelogram form 

T ■- \l,Et,E^t,...,E’^-h\ i (3.8) 

Let To be the leading n x n submatrix of T, and let Ti be the trailing n X n submatrix of T. Then 
To is a left triangular Toeplitz matrix, and T\ is a unit right triangular Toeplitz matrix. 



Example 3.1. Let n = 3. Then 



Xo 






Xi 


Xo 




X2 


X2 


Xo 


X3 


X2 


Xi 




X3 


X2 






X3- 



where we note that xs = L • 



To = 


Xo 

Xi 


Xo 


, Ti = 


X3 X2 
X3 


Xl 

X2 




_X2 


Xi Xo. 






X3. 



Lemma 3.2. Let To and T\ be defined as above. Then 

To" To + T"Ti = ToTo" + Ti T" . (3.9) 

Proof. Let N := T^T = TqTq + T/^Tj. We first show that is a Toeplitz matrix. Let Cj be 
defined by (3.1). Then by (3.8) we have for 1 < j\k < 

ej Nek = CjT^T E’^-H, 

where we have used that E^ = E~^ . We next define the reversal matrix 



J := [e„,e„_i,...,eij e 

Toeplitz matrices are counter symmetric, i.e. N = J J. Using that N^ is counter symmetric and 
Hermitian yields 

Tq"To + T"Ti = N = JN'^J = JNJ = J{T^Th + T[Ti)J 

= JTJ j ■ J%J + JT[ j • JT\J = ToTo" + TjT" . « 

The next lemma presents a Gaussian factorization of in terms of Tq and T\. This will be used 
together with Lemma 3.1 to express in terms of Tq and T\. 



Lemma 3.3. 






-ToTf'. 



(3.10) 



Proof. We first show that 






Vo 



VoA^ 



= 0 . 



(3.11) 



Let ey be defined by (3.7) and assume for the moment that the Xk are distinct. Then 

Vo 






VoA" 



Cfc = x{><k)X 



j-i 



(3.12) 
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and the right hand side vanishes for 1 < j\k < n. If the are confluent, then the right hand side 
expression of (3.12) contains derivatives of x('^) evaluated at A^. The right hand side of (3.12), 
however, still vanishes and (3.11) holds. 

We now write (3.11) as 

Tj'Vo + TfVoA" = 0 

and apply (3.2). This shows (3.10). • 

We are now in a position to show (1.18). By (3.4) G satisfies 

G = / + F"" 

and an application of the Sherman-Morrison- Woodbury formula yields 

Q-i = (/ + friQ ^ j _ F^{G~^ + , ( 3 . 13 ) 



We now determine an expression for 

r := I-G~\ (3.14) 

Substitute Y and (3.10) into (3.13) to obtain 



Y = To{T^To + Tl^Ti - T(^YTi)-^T^. (3.15) 

In order to determine a simple expression for Y from (3.15) we need the following observation, 
which is also central to Section 4. Tq and are both left triangular n X n Toeplitz matrices. 

Multiplication of Tq with can be identified with polynomial multiplication, see [Hel, Section 

1.3] and Section 4. Since multiplication of polynomials commutes, we obtain 

ToTi“^ = Tf^To. (3.16) 



From the correspondence between polynomials and left triangular Toeplitz matrices it also follows 
that TqT^^ is a left triangular Toeplitz matrix. 

Lemma 3.4. Equation (3.15) has the unique solution 

y = ri"^ToTo^Tf ^ = ToT-^Tf ^To^. (3.17) 

Proof. Unicity follows from (3.14) and that (3.4) hats a unique solution. From (3.16) we obtain 

= ToTr^Tr^T^- (3.18) 



Now substitute 
into (3.15). We obtain 



Y 






Tf "ToTo"Tf ' = To (To" To + T"Ti - To To")" 'To". 



(3.19) 
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An application of (3.9) reduces (3.19) to (3.18). The latter has already been shown to be valid. 
Therefore (3.17) solves (3.15). ■ 



Let 

Then Bq is a left triangular n x 
From (3.3) it now follows that 



Do := ToTj-^ = Tr^To. 
n Toeplitz matrix. By (3,14) and (3.17) 

G-' = I-BoBl = 1- BjBo. 
{R^R)-^ = I-BoB”. 



(3.20) 



(3.21) 
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4. Computation of Hq and Bq 

We summarize some results in [He 1, Section 1.3] and [He 2, Section 13.9] in order to show that the 
elements of Hq and Bq can be computed in 0{n log n) arithmetic operations from the coefficients 
Xj of X 3-nd 7Ty of 7T, see (l.l). To a polynomial or power series 

n-l 

sw ■■= E » 

J=0 

we associate the left triangular n x n Toeplitz matrix 

^ fy = 0 for j < 0, 

and we write f — > Z. If ^(A) is a polynomial and X a left triangular n X n Toeplitz matrix such 
that ^ > X, then it is easily seen that — > ZX. In particular, ZX is a left triangular n X n 

Toeplitz matrix. From and — > XZ is follows that ZX = XZ. 

Assume that f o 0 and let l/f — ► Z'. Then l/f • f — > /, Z'Z and ZZ'. We obtain Z' == Z“^ and 
therefore Z"^ is a left triangular Toeplitz matrix. 

Example 4.1. We have x ^o- Let 



X(A) ;= A" x(l/>) = Exn-,-''- (41) 

J=0 

Then x and the Blaschke product 

4 -ToTf" = 50. (4.2) 

X ■ 

Now let ^(A) and f(A) be arbitrary polynomials such that f (0) ^ 0. Henrici [He2, Theorem 13. 9e ] 
shows that the first n coefficients in the MacLaurin expansion of C(A)/f(A) can be computed in 
0(n log n) multiplications. The proof uses FFT. It is easily seen that the number of additions also 
is 0(n log n). 



From Xn = 1 and (4.1) we obtain x(0) ^ 0. Hence, the first n terms in the Ma cLaurin expansion 
of x/x can be computed in 0(n log n) arithmetic operations. By (4.2) therefore TqT^^ — Bq can 
be computed in 0(n log n) arithmetic operations. 



Because A^x(l/A) ^ 0 for A = 0, we can compute the first n terms in the MacLaurin expansion of 






in 0(n log n) arithmetic operations. This shows that Hq can be computed in 0(n log n) arithmetic 
operations. 
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5. A Factorization of R 



It follows from (3.3) and (3.20) - (3.21) that 

= Ir^~^ = I - Do = /-rr^ToTo^^T*\ 



and therefore 



T/^G”^Ti = T(^ 7 \ - ToT^ =: Afo"'. 



H 



(4.1) 

(4.2) 



The expression defining is a Gohberg-Semencul formula for the inverse of an n X n Toeplitz 

matrix, see, e.g., [lo, Theorem 18.2, p. 152]. We denote this Toeplitz matrix by Mq. From the left 
hand expression of (4.2) and the nonsingularity of Tj and R it follows that Mq is Ilerrnitian and 
positive definite. The desired factorization of R^ R is 



R^ R^T^ Mo 

We will now show how A/q can be computed. The computation involves running the Levinson 
algorithm backwards. 

Consider the related Gohberg-Semencul formula, see, e.g., (lo. Theorem 18.1, p. 148] or [AG], 

// 

Xn Xri— 1 • • • XO 



Mr^ = 



Xn-l 

Xn 



Xn Xn—1 



■ 0 




■ 0 


XO 




Xo 






Xl 


; 






-Xn-l ■■■ Xl Xo 0. 




-Xn-l ••• Xl Xo 0- 



XO 



Xn ~ 1 
Xn 






(4.3) 



where the four triangular Toeplitz matrices define the inverse of an (n + 1) X (n + l) Hermitian 
Toeplitz matrix. Denote this Toeplitz matrix by Mi. Then Mo is the leading principal n X n 
submatrix of Mi, see [lo. Theorems 18.1 - 18.2). 

Let Ri := \Pjk]'jk=zo ^ be the unit right triangular matrix, and let 

Di := diag[5o,Si, be the diagonal matrix such that 

= (4.4) 
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Given Mi — [f^j~k]^^k=oy matrices Ri and Di can be computed by the Levinson algorithm, and 
by comparing i?i with (4.3) one finds that 

Pjrt - Xjy 0 < y < n and (5^ = Xn, 

see, e.g., [AG]. We now apply the recursion formula in Levinson’s algorithm backwards in order 
to determine Ri and Di from the last column of Ri and (5,^. Then the recursion formula is used 
forwards to determine Mq. We will also obtain a Hermitian factorization of R^ R into triangular 
matrices. 



Backward Levinson algorithm 

input: \pjn\j^Qy^n ; output: Schur parameters {7y}y_i of Mq; 

for A: := n, n — 1, n — 2, ..., 1 do 

Ik := Pok\ Pk-i,k-\ := 1; 

for y := 1, 2, ..., integer part(|) do 

solve for and the linear system of equations 



■ 1 


7fc 




1 

»-» 

V 

1 




' _pj,k 




1 








. Pk-j,k . 



~ (&/(l - h»|))/(i + h»|); 

Levinson recursion for computing Mq = 
input: output: 

^lQ := So; /ii := -Sq^^; 

for k ~ l,2,...,n — 1 do 

Pk-i-l ~ ^^klfk-j-1 ~ X)y=l 



Hence Mq, i?i, and Di are computed in O(n^) arithmetic operations from the coefficients of x- Let 
Ro and Dq denote the n X n leading principal submatrices of i?i and Di respectively. Similarly to 
(4.4) we have 

R^MoRo = Do. (4.5) 

1 /2 

Because Mq is positive definite, so is Dq. Dq can therefore easily be computed. We obtain from 
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(4.1) - (4.2) and (4.5), with R := dI''^ 



R^R = { RTI^ f[RTl^). (4 .6) 

The right hand side of (4.6) is a Hermitian factorization into triangular matrices. It can be computed 
in 0{ n^) arithmetic operations from the coefficients of x- 
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